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The  results  of  a  packaqe  of  FORTRAN  computer 
programs  for  modelinq  defects  in  ionic  crystals  and  for 
fitting  experimental  data  are  described.  The 
fundamental  concept  of  the  defect  simulation  is  similar 
to  HADES  except  that  the  minimization  procedure  is 
different  since  the  package  is  designed  to  run  on  small 
computers.  As  an  example  of  the  use  of  this  packaqe, 
the  relative  stabilities  of  nn  and  nnn  complexes  for 
various  rare  earths,  lanthanum,  and  yttrium  are 
considered.  First,  the  data  fittinq  routine  was  used 
to  analyze  relaxation  data  for  nn  and  nnn  complexes  in 
rare  earth  doped  strontium  fluoride.  The  experimental 
results  for  strontium  fluoride  were  then  used  in 
conjunction  with  che  defect  simulation  program  to 
determine  potentials  for  all  of  the  rare  earths, 
yttrium,  and  lanthanum.  Those  rare  earth  potentials 
were  then  used  in  the  simulation  of  calcium  fluoride 
and  shiw  that  the  nnn  complex  should  not  be  observable 
except  for  possibly  the  smallest  rare  earths.  This 
implies  that  the  B  site  of  Wright  and  co-workers  or  the 
RII  relaxation  requires  another  explanation.  Also,  the 
potentials  were  used  in  the  simulation  of  barium 
Muoride,  showinq  that  the  nn  complex  should  be 
observable  only  for  the  largest  rare  earths  or 
lanthanum.  Next,  the  enthalpy  for  nn-»nn  reorientation 


via  the  i nter st i t i a  Icy  mechanism  was  calculated  for 
rare  earths  in  calcium  and  strontium  fluoride.  In 
general,  the  calculated  reorientation  enthalpies  are 
larger  than  the  experimental  values.  However,  the 
variation  of  the  enthalpy  with  the  size  of  the  rare 
earth  is  in  reasonable  agreement  with  experiment. 
Finally,  the  variation  of  the  calculated  enthalpy  with 
pressure  is  found  to  be  in  excellent  agreement  with 
exper iment . 


1.  Introduction 


In  the  past  few  years,  the  authors  have  been 
studying  defects  in  ionic  crystals  (Andeen  et  al,  1981; 
Andeen  et  al,  1977),  using  dielectric  relaxation 
techniques.  The  relaxation  spectra  of  many  systems 
previously  thought  to  be  quite  simple,  are  found 
actually  to  be  very  complex.  For  example,  seven 
relaxations  are  observed  in  rare  earth  doped  calcium 
fluoride  (Andeen  et  al.  1981). 

In  order  to  aid  in  the  interpretation  of  the 
data,  the  authors  have  developed  computer  simulation 
techniques  similar  to  the  HADES  (Harwell  Automated 
Defect  Evaluation  System)  program,  which  has  produced 
many  significant  results  (Catlow,  Norgett,  and  Ross 
1977;  Catlow  1976).  The  main  concept  is  to  model  a 
crystal  by  a  central  region  of  movable  ions  surrounded 
by  a  dielectric  continuum.  Shell  model  potentials  are 
used  for  the  movable  ions.  Defects  may  then  be 
introduced  into  the  central  region  after  which  the  ions 
are  allowed  to  relax  to  their  minimum  energy 
configuration.  However,  the  difficulty  with  HADES  is 
that  the  storage  requirements  are  great,  requiring  a 
very  large  and  fast  computer.  Consequently,  the 
authors  have  developed  an  alternative  package  which 
requires  far  less  storage  with  only  a  modest  increase 
in  the  running  time,  making  it  suitable  for  more 
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general  use. 

The  particular  system  used  by  the  authors  has 
the  additional  feature  that  all  output  may  be  displayed 
on  an  Evans  and  Sutherland  Picture  System. 

Stereographic  display  is  available  to  enhance  the 
three-dimensional  effect.  Movie-making  facilities  also 
tie  in  with  the  picture  display. 

As  a  related  application  of  the  minimization 
procedure  associated  with  the  computer  simulation,  a 
program  for  the  least-squares  fitting  of  data  has  also 
been  developed.  In  the  example  presented  in  this 
paper,  both  techniques  are  used  to  study  simple  point 
defects  in  rare  earth  doped  alkaline  earth  fluorides. 

2 •  The  Programming  Technique 

The  key  feature  in  computer  programs  of  this 
type  is,  of  course,  the  minimization  procedure.  In  the 
case  of  data  fitting,  the  sum  of  the  squares  of  the 
differences  between  the  data  and  the  best-fit  curve  is 
minimized.  In  simulating  a  defective  solid,  the  energy 
of  the  central  region  will  be  minimized.  The  following 
brief  discussion  of  the  minimization  procedure  will  be 
in  terms  of  the  defect  simulation.  A  similar  technique 
is  used  in  data  fitting,  where  adjustments  are  made  in 
the  fitting  parameters  rather  than  ion  positions. 

In  the  present  program,  a  modified  block  guasi- 
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Newton  method  (Rheinboldt  1974)  is  used.  By  block,  we 
mean  an  iterative  approach,  where  the  position  of  each 
ion  is  adjusted,  one  at  a  time.  This  is  the  analogue 
of  the  Gauss-Seidel  procedure  for  solving  a  linear 
system.  The  position  of  each  ion  is  represented  by  six 
variables,  three  nuclear  and  three  shell  coordinates. 

A  standard  Newton-Raphson  method  is  used  to  determine  a 
new  location  for  the  chosen  ion,  except  that  for 
stability  of  the  algorithm,  and  as  a  safeguard  against 
local  minima,  the  energy  is  evaluated  at  eight 
positions  nearby.  These  values  are  used  to  approximate 
the  derivatives  called  for  in  the  Newton-Raphson 
technique.  There  are,  of  course,  risks  involved  in 
numerical  approximation  of  derivatives  and  great  care 
is  taken  to  avoid  the  pitfalls. 

However,  convergence  using  this  procedure  is 
slow.  Thus,  the  procedure  was  modified  to  evaluate  the 
potential  function  twice  in  the  direction  of  the 
Newton-Raphson  step,  at  the  Newton  distance  and  at  half 
that.  Fitting  to  a  parabola  determines  much  more 
accurately  the  minimum  in  this  direction.  As  a  check, 
the  energy  is  evaluated  once  more  at  the  location 
determined  by  the  parabolic  fitting.  The  ion  is  moved 
to  the  position  at  which  the  energy  was  least.  This 
provides  a  safeguard  against  error  in  the  numerical 
approximation  of  derivatives.  R-linear  convergence 


(Rheinboldt  1974)  is  expected  and  this  has  been 
observed . 

Symmetry  is  not  invoked  in  our  minimization 
and  thus  configurations  of  low  symmetry  are  easily 
dealt  with  by  the  current  program. 

The  most  important  feature  of  the  present 
implementation  is  that  the  program  requires  only  18,500 
3fi-bit  words  together  with  the  ability  to  perform 
double  precision  computations  in  FORTRAN.  This, 
combined  with  the  efficient  minimization  procedure 
makes  the  program  suited  to  small  computers. 

In  the  case  of  data  fitting,  at  first  inspection 
it  might  appear  that  a  more  efficient  minimization 
technique  such  as  that  of  Broydon  or  Fletcher-Powell 
would  be  more  desirable  since  the  storage  requirements 
would  be  small.  However,  the  present  technique  does 
not  require  calculation  of  the  derivative  of  the 
function  to  be  fit.  In  addition,  the  difference  in 
times  involved  in  minimization  by  the  different 
techniques  is  small  because  of  the  the  number  of 
variables  involved.  Finally,  the  fits  are  done  in  an 
iterative  mariner  and  the  resultant  graphical  display 
provides  added  evidence  of  the  goodness  of  fit. 
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Consequently,  the  minimization  technique  employed  in 
the  defect  simulation  was  also  applied  to  the  data 
fitting . 

Further  details  of  the  techniques  are  given  in 
Appendices  A-C. 

3 .  Defect  Simulations 

We  now  give  further  details  of  the  operation  of 
the  defect  simulation  package.  The  model  contains 
about  150  movable  ions  (5x5x5  fluorines  plus 
appropriate  cations) .  These  are  "surrounded"  by  a 
continuum  of  fixed  ions,  each  of  which  is  assigned  a 
total  (electronic  plus  ionic)  polarizability.  The 
electrostatic  calculations  for  a  given  ion  in  the 
central  region  are  performed  exactly  for  the 
surrounding  cube  of  llxllxll  fluorines  plus  cations. 
The  remaining  electrostatic  energy  is  approximated  via 
a  Madelung  sum. 

One  or  more  defect  ions  can  be  added  to  a 
computer  file  containing  the  basic  model  including 
subst i t ut iona I s ,  i nterst i t ial s ,  or  vacancies.  The 
lattice  is  then  allowed  to  relax  to  the  minimum  energy 
configuration  as  described  in  the  previous  section. 

The  energy  of  a  given  ion  is  the  sum  of  five 
terms,  E=RSP+FS+FNtFNR+ENR2 .  The  spring  energy,  FSP, 
is  due  to  the  relative  positions  of  the  nucleus  and 


shell.  The  shell  component  of  the  electrostatic  enerqy 
is  ES  and  EN  is  the  nuclear  component.  The  near 
neighbour  repulsive  interaction,  ENR,  is  given  by 
summing  the  potentials  A  exp(-r/  P)  over  the  nearest 
neighbours  of  the  ion  being  moved.  A  and  P  are 
constants  depending  on  the  ions  involved.  The  second 
neighbour  interactions,  ENR 2 ,  are  included  in  exactly 
the  same  manner  as  by  Catlow  et  al.  (1977). 

The  distance  between  ions  is  limited  to  be  no 
smaller  than  1.25  times  one  lattice  unit  (defined  as 
half  the  average  f luor ine-f luor ine  distance),  so  as  to 
prevent  computational  instability  in  the  form  of  two 
ions  "collapsing  together."  In  each  case  the  program 
searches  for  the  two  ions  having  the  smallest 
separation  after  final  convergence.  In  no  case  have 
the  two  closest  ions  been  near  1.25  lattice  units 
apart,  implying  that  the  arbitrary  repulsion  has  not 
affected  the  results. 

In  the  present  model,  the  continuum  polarization 
was  not  included  in  the  minimization  procedure. 
Certainly,  this  has  little  effect  in  the  case  of 
neutral  defects.  However,  in  order  to  evaluate  the 
magnitude  of  this  effect,  the  polarization  enerqy  was 
calculated  after  minimization,  as  follows. 

The  excess  charge,  0,  and  dipole  moment,  P,  of 
the  central  movable  region  was  evaluated.  Then  the 


energy  of  each  ion 
to : 


was  set  equal 


a.  (e'-l) 

w  - - i - - - 

a  (4m:  ’  )  iQNr 


(Q2r2+4PQr (cosO) +P2+3P2cos2 0 } 


where  e '  is  the  static  dielectric  constant  and  N  is 
the  number  of  molecules  per  unit  volume.  This  is 
essentially  a  modified  Mott-Littleton  (1938)  approach 
since  the  polarizabilities,  defined  by: 


and 


e'+2  3c 

o 

+  2ot,., 

anion  fluorine 


(2) 

(3) 


include  both  electronic  and  ionic  polarizabilities. 
Further,  the  polarizabilities  were  separated  usinq  the 
empirical  result  that: 


when1  >  is  the  isothermal  compress!  bli  ty . 

This  result  is  known  as  the  Jarman  rule  (Gibbs  and 
orman  19^2)  and  i he  experimental  results  for  the 
alkaline  earth  fluorides  usinq  the  data  of  Andeen  et 
al.  (1972)  are  listed  in  table  1.  Also  listed  there 
are  the  values  for  some  other  ionic  crystals  usinq  the 
data  of  Fontanel  la  et  al.  (  1972).  ft  is  seen  that  the 


above  relation  holds  quite  well.  Since  equation  4 
implies  that: 

2 

a  =  (Constant)V 

it  is  assumed  that: 

a  =3  (Constant) 

Using  the  ionic  radii  of  Shannon  (197^),  a  set  of 
polarizabilities  was  calculated  and  these  are  listed  in 
table  2. 

4 .  Simple  Point  Defects  in  Rare  Earth  Doped  Alkaline 

Earth  Fluorides 

4.1.  Data  Fitting 

Simple  point  defects  in  rare  earth  doped 
strontium  fluoride  provide  an  excellent  case  in  which 
to  apply  data  fitting  and  defect  simulation  techniques. 
This  can  be  seen  in  figure  2  of  a  previous  paper 
(Andeen  et  al.  1981)  where  it  is  seen  that  for 
strontium  fluoride  the  ratio  of  the  two  principal  peaks 
varies  smoothly  with  a  change  in  the  size  of  the  rare 
earth.  What  is  not  apparent  is  that  in  cases  where 
both  peaks  exist  simultaneously,  they  are  in  thermal 
equilibrium  with  one  another.  This  is  clearly  shown  in 
fiqure  1  where  the  results  for  gadolinium  doped 
strontium  fluoride  are  plotted.  It  is  seen  that  the 
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height  of  the  higher  temperature  peak  actually 
increases  with  temperature  while  the  height  of  the 
lower  temperature  peak  decreases  faster  than  expected 
on  the  basis  of  Curie-Weiss  behaviour.  This  implies 
that  the  higher  temperature  peak  grows  at  the  expense 
of  the  lower  temperature  peak.  This  problem  has,  of 
course,  been  solved  by  Matthews  and  Crawford  (1977)  who 
explain  their  ITC  data  for  gadolinium  doped  strontium 
fluoride  by  assigning  the  peaks  to  jumps  of  nn  and  nnn 
dipoles.  They  assign  an  enthalpy  difference: 


AE 


12 


E 


12 


E 


21 


to  the  difference  in  binding  energies  of  the  two 
complexes.  This  corresponds  to  the  enthalpies  shown  in 
figure  2,  and  AE12  will  be  a  key  parameter  in  the 
following  discussion. 

The  model  of  Matthews  and  Crawford  was  applied 
to  the  data  in  the  following  manner.  Firstly,  Cole- 
Cole  expressions  were  used  to  approximate  the 
relaxation  peaks  leading  to: 

2  A  .  cos  (a',  n/2 ) 

..  _  v  _ }_ _ i_ _ 

■  2T!cosh(  ( 1— nt  t  )  x .  )  +sin  (<x  ti/2  )  I 

l— I  li  l  (51 


where 

x ^  In  (i.i  (  io)ioxp(Eil/kT)  ) 
12 


T — 


a ! 

1  is  the  Cole-Cole  parameter  and  is  a  measure  of  the 

]•;  jr 

broadeninq  of  the  peaks.  The  use  of  11  and  ‘21 
follows  from  the  approximation 


J11 


E12  and  E21 


(  t  )  . 

as  shown  by  Matthews  and  Crawford  (1977).  The  °  1 
are  the  so-called  reciprocal  frequency  factors.  Also, 


A  =  -1— X- 
i  3e  k 
o 

N.  p 

where  i  is  the  dipole  concentration  and  i  is  the 

dipole  moment.  It  follows  that: 


exp(AS12/k  -  AE12/kT) 


(7) 


Conservinq  the  total  number  of  dipoles: 
N1  +  N2  =  NT 

therefore : 


I  .1 +oxp  ( AS  ^  2/k  -  AH  |  2/kT)  I 


and 


A, 


(  l  +  oxp(AE12/kT 


AS12/k)  1 


(9) 


n 


Consequently,  the  theoretical  expression  given  by  the 
above  equations  contains  ten  adjustable  parameters. 


AS  12 ,  AE 


12 


This  is,  of  course,  a  large  number  of  parameters. 
However,  it  should  be  remembered  that  effectively  ten 
peaks  are  being  fitted  simultaneously  since  there  are 
two  different  relaxations  with  five  curves 
(frequencies)  for  each.  (Only  three  are  shown  in 
figure  1.) 


The  results  of  the  two  peak  data  analysis  are 
given  in  table  3  for  the  central  four  rare  earths.  A 

typical  best  fit  curve  is  shown  with  the  data  in  fiqure 

AE 

1.  The  value  of  12  for  gadolinium,  0.058  eV,  is 
slightly  larger  than  the  values  0.045,  0.049,  and  0.044 
eV  reported  by  Matthews  and  Crawford  (1977),  Fontanella 
et  al.  (1978a),  and  Edgar  and  Welsh  (1979), 
respectively.  However,  they  probably  agree  to  within  a 
reasonable  error.  The  value  of  0.00  eV  for  europium  is 
close  to  0.078  eV  reported  by  Fontanella  et  al. 

(1978a).  The  other  values  are  new.  The  peaks  for  the 
other  materials  were  determined  using  single  peak  fits 
to  ocpi.it ion  (')  )  . 
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4.2.  Computer  Simulations 
4.2.1.  Defect  Stabilities 

The  computer  simulation  was  then  applied  to  the 
same  problem.  The  interionic  potentials  are  the  same 
as  those  used  by  Catlow  et  al.  (1977).  As  an  initial 
test  of  the  model,  a  suhstitutional  trivalent  rare 
earth  was  modelled  merely  by  adjusting  the  shell  charge 
of  the  calcium  ion,  in  the  same  way  that  Catlow  (1975) 
did.  In  this  approximation,  it  was  found  that 

(AE12)th  =  -0.11  and  -0.31  eV 

for  calcium  and  strontium  fluoride,  respectively. 

These  values  are  in  reasonable  agreement  with  the 
values  of 

(AE12)th  =  -0.03  and  -0.4  eV 

determined  by  Catlow  (1975)  using  the  HADES  model.  In 
part,  the  small  difference  is  due  to  the  polarization 
of  the  continuum.  However,  the  important  result  is 
that  the  theoretical  values  are  much  lower  than  the 
experimental  values  listed  in  table  3  since  the  values 
for  Tb-Sm  are  all  positive.  Since  the  value  of  ^^12^EXP 
decreases  as  the  size  of  the  rare  earth  decreases,  it 
was  concluded  that  the  error  in  the  theoretical 
calculation  was  with  the  rare  earth  potential,  the 


effective  ionic  radius  being  too  small. 


The  simplest  approach  was  taken  to  simulate 
larger  trivalents,  namely,  the  pre-exponential  in  the 
expression  for  the  near-neiqhbour  interaction 

V  (r)  =  A  exp  (-r/p) 

was  increased.  The  results  of  the  calculations  for 
relative  binding  enerqies  using  various  values  of  A  are 
listed  in  table  4.  The  best  fit  straight  line  for 

^AE12)TH  VS  A  WaS  determined*  Using  (Ar]2*exp'  ifc  was 
concluded  that  Sm,  Eu,  Gd,  and  Tb  correspond  to  values 
of  A  of  2488,  2408,  2349,  and  2251  eV,  respectively. 

The  average  interval  is  about  79  eV  per  rare  earth. 
Extrapolating,  then,  values  for  A  for  all  the  rare 
earths  were  determined  and  listed  in  table  S.  Also, 
yttrium  is  about  the  same  size  as  erbium  (Fontanella  et 
al.  1978b;  Shannon  1978)  so  that  A=?014  eV. 

One  of  the  advantages  of  the  interface  with  the 
3D  picture  system  is  that  it  is  possible  to  see  the 
reason  for  the  variation  of  the  relative  stability  of 
the  two  complexes.  This  is  shown  in  figures  3  and  4 
where  orthographic  projection  plots  of  nn  and  nnn 
dipoles  are  shown  for  an  intermediate  (A=2300  eV)  and  a 
small  (A=18(10  eV)  size  trivalent  ion  in  strontium 
fluoride.  Clearly,  the  substitutional-interstitial 
interionic  distance  is  smaller  for  the  smaller 
trivalent  for  the  nnn  complex  but  is  actually  larger 


for  the  smaller  trivalent  for  the  nn  complex.  This  is 
an  example  of  distortion  causinq  nn  complexes  to  become 
less  stable  relative  to  the  nnn  complex  as  the  size  of 
the  rare  earth  decreases.  That  A=1273  eV  corresponds 
to  a  very  small  trivalent  ion  can  be  seen  in  figure  5  . 
Other  examples  of  the  output  of  the  3D  picture  system 
are  given  in  figures  3  and  4  of  the  paper  by 
Wintersgill  et  al.  (1980)  and  fiqures  3  and  4  of  the 
paper  by  Andeen  et  al.  (1981). 

Following  the  determination  of  potentials  for 
various  rare  earths 

in  strontium  fluoride  those  potentials  were 

used  in  simulations  of  calcium  fluoride  and  barium 
fluoride.  The  results  of  the  calculations  are  listed 
in  table  4. 

The  results  for  calcium  fluoride  can  be  used  to 
comment  on  a  recent  controversy  concerning  the  RII 
relaxation  or  the  R-site  of  Wright  and  co-workers 
(Moore  and  Wright  1979,  1981;  Tallant  and  Wright  1975; 
Tallant  et  al.  1977).  It  has  been  shown  (Fontanella  et 
al.  1980)  that  the  relaxation  data  for  the  RII  peak 
lead  to  difficulties  concerning  the  relaxation 
mechanism  if  it  is  assigned  to  a  nnn  complex. 
Specifically,  equilibrium  studies  indicate  that  RII 
cannot  be  associated  with  a  nnn-fnn  jump.  That  is 
because  the  low  activation  energy  (0.15  eV)  implies 
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that  the  nn  and  nnn  sites  should  equilibrate  quickly  at 
room  temperature  and  thus  remain  in  the  same  ratio. 

•The  experimental  data  do  not  support  this  conclusion. 
The  only  alternative  is  a  nnn->nnn  reorientation 
mechanism;  however,  the  very  low  activation  energy 
(0.15  eV)  argues  against  this  conclusion. 

The  results  of  the  present  work  support  the 
assertion  that  RII  or  the  B-site  cannot  be  assigned  to 
a  nnn  complex.  That  is  because  the  calculations  for  a 
rare  earth  the  size  of  erbium  (A=2014  ev;  yield 

(AE12)TH  =  +0'21  eV 

Consequently,  at  100K,  about  the  position  of  RII,  only 
about  2  in  10 ^  of  the  ions  associated  with  the  simple 
point  defect  will  be  in  the  nnn  position.  This 
conclusion  is  supported  by  Baker  (1974)  who,  in  his 
review  of  the  ESR  literature  implies  that  the  nnn 
complex  has  never  actually  been  observed  in  rare  earth 
doped  calcium  fluoride  using  ESR  techniques. 

It  is  interesting  that  RII  is  only  observed  for 
small  rare  earths  (Andeen  et  al.  1977).  This  is,  of 
course,  the  behaviour  which  might  be  expected  for  a  nnn 
complex  on  the  basis  of  the  present  calculations  though 
the  magnitude  of  the  site  populations  are  not 
consistent  with  this  identification.  However,  this 
same  behaviour  is  observed  for  RII I  which  has  been 
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correlated  with  the  D(2a)  site  of  Wriqht  and  co-workers 
(Fontanella  et  al.  1980;  Andeen  et  al.  1981)  and  which 
is  some  sort  of  cluster.  In  fact,  it  has  been 
suggested  (Andeen  et  al.  1981)  that  this  intimate 
relationship  between  RII  and  RIII  may  lead  ultimately 
to  the  identification  of  these  two  complexes. 

The  simulation  results  for  barium  fluoride  are 
listed  in  table  4,  and  predict  that  the  nn  complex  will 
only  be  stable  relative  to  the  nnn  complex  for  the 
largest  rare  earths  or  lanthanum.  This  is  observed 
exoer imental ly  (Laredo  et  al.  1979,  1980;,  Andeen  et 
al.  1981),  a  low  temperature  relaxation  region  being 
present  for  lanthanum  and  cerium  (nn)  and  and  a  higher 
temperature  relaxation  region  (nnn)  beinq  observed  for 
all  of  the  rare  earths.  However,  the  higher 
temperature  relaxation  region  consists  of  two  closely 
spaced  peaks  (Laredo  et  al.  1980;  Andeen  et  al.  1981). 
Laredo  et  al.  (1980)  have  identified  one  of  those  as  a 
nn  relaxation  which  they  observe  for  small  rare  earths. 
Clearly,  the  present  calculations  suggest  that  this  is 
unlikely  as  there  is  no  tendency  for  the  nn  complex  to 
become  stable  relative  to  the  nnn  complex  for  small 
rare  earths. 

As  pointed  out  previously  (Andeen  et  al.  1981), 
the  conclusion  that  only  nnn  complexes  exist  in  barium 
fluoride  doped  with  small  rare  earths  is  supported  by 
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the  selective  laser  excitation  work  of  Miller  and 
Wright  (1978a,b)  in  that  they  find  only  one  type  of 
simple  site  in  erbium  doped  barium  fluoride. 

Clearly,  then,  a  nn-»nn  reorientation  mechanism 
cannot  be  responsible  for  one  of  the  peaks  in  barium 
fluoride  doped  with  small  rare  earths.  It  is  tempting 
to  attribute  the  two  peaks  to  nnn— »nnn  and  nnn-tnn 
reorientations.  However,  that  cannot  he  the 
explanation  as  both  peaks  would  not  be  observed  via 
ITC. 

A  better  explanation  for  one  of  the  two  peaks  in 
barium  fluoride  doped  with  small  rare  earths  is  a 
cluster.  In  particular,  the  position  of  one  of  the 
peaks  does  depend  moderately  on  the  size  of  the  rare 
earth  (Laredo  et  al.  1980).  In  fact,  the  RTV 
relaxation,  which  has  been  observed  in  both  calcium  and 
strontium  fluorides  (Andeen  et  al.  1977,  1981)  and 
which  has  been  associated  with  a  "gettered" 

(Wintersgill  et  al.  1980;  Andeen  et  al.  1981),  depends 
strongly  upon  the  nature  of  the  rare  earth.  However, 
the  experimental  data  concerning  the  barium  fluoride 
relaxations  to  date  (Laredo  et  al.  1980;  Andeen  et  al. 
1981)  do  not  support  assigning  either  to  a  cluster. 
Further  work  is  necessary  for  their  identification. 


.  Effect  of  Pressure  on  Reorientation  Enthalpies 


As  a  second  test  of  the  simulation  proqram,  the 
reorientation  enthalpy  was  calculated  for  jumps  of  an 
interstitial  between  nearest  neighbour  equivalent  sites 
via  the  interstitialcy  mechanism  for  calcium  fluoride 
and  strontium  fluoride.  In  the  present  calculations, 
the  saddle  point  of  the  migratinq  ion  was  not 
predetermined.  Only  one  coordinate  of  the  diffusing 
ion  was  fixed.  The  saddle  point  was  then  determined  by 
minimization  at  successive  values  of  a  single 
coordinate  of  the  diffusing  ion.  The  activation 
enthalpy  was  then  taken  to  be  the  maximum  value  of  the 
lattice  energy  along  the  path  of  the  diffusing  ion. 

The  results  are  listed  in  table  6  together  with  the 
corresponding  experimental  values.  The  theoretical 
values  are  about  0. 1-0.2  eV  higher  than  the 
experimental  values.  While  it  is  temptinq  to  assign 
the  discrepancy  to  the  rare  earth  potential,  the 
difference  may  be  more  fundamental  in  that  there  is  no 
reason  to  expect  a  static  calculation  to  accurately 
simulate  an  inherently  dynamical  process.  However,  the 
calculations  do  show  only  a  small  decrease  in  the 
enthalpy  as  the  size  of  the  rare  earth  decreases.  That 
is  in  agreement  with  experiment. 

Next,  the  effect  of  pressure  on  the 
reorientation  enthalpy  for  nn-Vnn  jumps  in  strontium 
fluoride  doped  with  a  large  rare  earth  was  determined 
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theoretically  from 


=  -XTV 


where  T  is  the  isothermal  compressibility  for 
strontium  fluoride.  The  enthalpy  was  found  to  increase 
by  about  0.0042  eV  when  the  lattice  constant  was 
decreased  by  about  0.05%  (0.1  GPa  pressure).  Thus, 


(f§  =42  meV/GPa 

'  d  'th 


which  is  in  excellent  agreement  with  the  experimental 
value  of  41  meV/GPa  (Fontanella  et  al.  1981).  (In  the 
present  calculation,  the  central  region  was  not  further 
relaxed  upon  application  of  pressure,  only  the  lattice 
constant  was  changed.)  The  agreement  is  better  than 
expected  considering  the  uncertainty  in  the  experiment 
and  the  nature  of  the  theoretical  calculation. 

In  summary,  while  the  program  does  not  yield 
good  values  of  reorientation  enthalpies,  the  effects  of 
small  perturbations,  such  as  pressure  or  small  shifts 
in  ion  size,  on  the  enthalpy  are  reproduced  quite  well. 


5 .  Conclusions 

In  summary,  then,  the  results  of  a  new  computer 
program  for  modeling  ionic  solids  and  best-fitting  data 
are  presented.  The  techniques  are  applied  to  the 


problem  of  simple  point  defects  in  rare  earth  doped 


alkaline  earth  fluorides.  A  combination  of  the 
techniques  has  led  to  the  following  results: 

1.  Potentials  are  obtained  for  all  the  rare  earths, 
lanthanum,  and  yttrium.  Those  were  obtained  by  a 
combination  of  theoretical  and  experiments]  results  for 
rare  earths  in  strontium  fluoride. 

2.  Calculations  for  rare  earths  in  calcium  fluoride 
clearly  show  that  the  nnn  complex  should  not  be 
observable  in  the  vicinity  of  room  temperature  or  below 
except  for  possibly  the  smallest  rare  earths  and  thus 
rule  out  assigning  the  RII  relaxation  or  the  B  site  to 
a  nnn  complex. 

3.  Calculations  for  barium  fluoride  clearly  show  that 
the  nn  complex  should  only  be  observable  for  the 
largest  rare  earths  or  lanthanum  and  thus  rule  out  nn-» 
nn  reorientations  for  small  rare  earths. 

4.  Calculated  nn— »nn  reorientation  enthalpies  for 
calcium  fluoride  and  strontium  fluoride  are  found  to  be 
larger  than  the  experimental  values,  however,  the 
variation  of  the  enthalpy  with  the  size  of  the  rare 
earth  and  with  pressure  agree  reasonably  well  with 

':•<  per  iment . 
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Table  1.  Values  of  the  logarithmic  derivative  of  the 
polarizability  determined  using  the  data  of 
Andeen  et  al.  (1972)  and  Fontanella  et  al.  (1972). 

V 

a 


CaF2 

2.09 

SrF2 

2.03 

BaF2 

1.92 

LiF 

2.01 

NaF 

2.32 

NaCl 

2.12 

NaBr 

2.05 

KC1 

2.05 

KBr 

2.00 
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**«*W  V 


Table  2 


Polarizabilities  of  ions  in  the  alkaline 

-40  2 

earth  fluorides  in  units  of  10  N/Coul  -m. 
The  values  include  both  electronic  and  ionic 
effects . 


aF 

aAE 

Ca 

1.96 

5.02 

Sr 

1.66 

6.54 

Ba 

1.39 

9.17 
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Table  4.  Results  of  the  computer  simulation  for  <AEi2^TH 


A  (eV) 

(AE12) 

™<eV) 

CaF2 

SrF2 

BaF2 

2700 

+  0.55 

+  0.19 

-0.022 

2600 

+  0.50 

+  0.16 

-0.057 

2400 

+  0.41 

+0.079 

-0.128 

2200 

+  0.31 

-0.002 

-0.20 

2000 

+  0.20 

-0.08 

-0.26 

1800 

+  0.09 

-0.15 

-0.29 
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Table  5.  Values  of  A  in  the  rare  earth  potential 


Table  6.  Reorientation  enthalpies  for  nn-*nn  jumps  of  an 
interstitial  via  the  interstitialcy  mechanism. 


CaF2 

A  (eV) 

Theory 

2600 

i^/Nd) 

0.69 

2400 

(~Eu) 

0.67 

2000 

( ***  Er ) 

0.61 

SrF2 

2600 

( Nd) 

0.564 

2400 

( r->  Eu) 

0.562 

a.  Andeen  et  al.  (1977). 

b.  From  table  3. 


i2 


Experiment 
0.4  3a 
0 . 414a 
0 . 406a 

0 . 4  8  b 
0.4  6b 


Figure  3.  Orthographic  projection  along  the  (100)  direction  for  an  intermediate  size  rare  earth  in 
strontium  fluoride  (A=2300  eV)  .  The  plus  indicates  the  position  of  the  rare 
earth  and  the  tetrahedron  represents  the  charge  compensating  interstitial 
fluorine.  The  straight  lines  connect  lattice  fluorines.  (a)  nn  complex 


Appendix  A.  The  Minimization  Philosophy. 

The  minimization  strategy  is  the  same  in  both 
our  parameter  fitting  of  data  and  in  the  computer 
modeling  programs.  While  the  basic  idea  is  well  known, 
its  implementation  in  this  application  is  not. 

In  order  to  develop  the  notation,  let  us  first 
consider  the  question  of  finding  a  solution  to  the 
vector  equation: 

g(a)  =  0 


where 


a  -  l‘Va2,--  'V' 


and 


g(a)  =  (g,  (a)  ,g_  (a)  ,  . .  .  ,g  (a)  J . 

^  ~  JL  ~  £  —  p 

Thus,  we  are  trying  to  find  a  simultaneous  solution  to 
the  p  nonlinear  equations: 


gi (ai'a2' • *  * ,ap*  =  °* 

Suppose  we  refer  to  such  a  solution  as 

a*  =  {af , a* , . . . ,a*  ) 

rsj  1  c.  P 

and  we  are  given  a  vector,  a,  which  is  close  to  our 

A/ 

solution,  a*.  Then  in  order  to  find  our  solution,  we 
have  just  to  find 

Aa  =  a*  -  a  or  Aa .  =  a?  -  a.  . 

If  Aa  is  sufficiently  small,  then  vector  calculus 
assures  us  of  the  approximat ion 

gj(a  +  Aa)  ~  gj(a)  ^ 

=  9 j <a)  +  (5Sj) • (Aa) « 

where 

Dg  .  Dg . 

V  g  .  =  f  -w — 1  r—  J  ) 

Da'  'Da 
J  1  P 

The  above  equations  can  be  summarized  by  the  equation: 
g(a+Aa)«  g(a)  +  (J  )  (Aa)  (1) 
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where  is  the  p  by  p  Jacobian  matrix  given  by: 


Of  course,  we  want  Aa  to  satisfy  g(a+  Aa)=o.  Thus,  we 


<v  <v 


use  equation  (1)  to  solve  for 


Aa  =  -  ( J  )  1 (g (a) ) . 


(2) 


The  technique  for  finding  a*  then  amounts  to 

a 

finding  a  close  initial  value,  say  ~°  ,  and  then 
iterating  by  the  formula 


a,,  =  a  -  ( J  )  ~ 1  (g  (a  )  ) 
^n+J.  »>»n  ^ Q 


(3) 


a  a 

where  the  Jacobian  is  evaluated  at  ~n  .  If  ~°  is 
"sufficiently  close,"  then  the  sequence 


~o'~l'  ~n 

* 

will  converge  fairly  rapidly  to  the  solution, 

This  is  of  course  the  Newton-Raphson  iterative  method 
for  finding  a  minimum.  In  the  computer  model,  we  let  S 
represent  the  total  energy  of  the  ions.  For  data 
fitting,  S  is  h:e  sum  of  the  squares  of  the  differences 
between  the  data  and  fitted  values. 


> ! >1  »■  miix  B.  Curve  Fitting. 

For  data  fitting,  wo  are  given  data  values: 


t  x  .  ,  y  .  1 

l  1  l 
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and  the  equation  to  be  fitted: 


y  =  f(x,£) 


where  a  represents  the  set  of  parameters  to  be  fit.  A 
solution  is  merely  a  value,  a*,  that  minimizes  the 
value  of  the  function: 

N  ? 

S  -  X  jy  -  f  (x  , n)  (4) 

i  =  l 

where  N  is  the  number  of  data  points.  In  this  case,  we 
are  attempting  to  find  a  solution  to  the  vector 
equation : 


VS 

<*»* 


r  3S  ,  3S  , 

3a^  ^^2 


,3S 

3a 

P 


0 


.  a 

Given  a  starting  value,  t  we  then  iterate  using  the 
appropriate  modification  of  equation  (3): 


:n+l 


=  a 
■~n 


-  (H)_1(VS) 


(5) 


where  li 

S ;  i  .  e . 


is  the  (less i an  (matrix  of  second  partials)  of 


H 


32S 


3a . 3a  . 
i  1 


and  is  the  gradient,  both  evaluated  at  ~'n  . 

For  our  curve-fitting,  we  have  developed 
software  using  APL  on  the  local  timesharing  system, 
NATS  (Naval  Academy  Time-Sharing  System),  which  is 
based  on  the  PTSS  (Partmouth  Time  Sharing  System).  We 
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chose  APL  for  the  following  reasons: 

1.  APL  is  an  interpreter  and  therefore  lends 
itself  well  to  interactive  programming. 

2.  APL  has  a  wide  range  of  built-in  functions 
including  virtually  all  of  the  so-called  scientific 
functions. 

3.  Our  implementation  of  APL  has  excellent 
built-in  graphics  capabilities. 

The  APL  package  we  have  developed  has  routines 

for  data  input,  least  squares  fitting  and  graphics 

display  of  the  data  along  with  the  theoretical  curve 

given  by  the  current  values  of  the  fitted  parameters. 

Some  details  about  the  least  squares  routine  are 

now  in  order.  One  of  the  functions  included  in  the 

package  is  a  function  S  that  depends  only  on  the 
a  -1  a 

parameters  x  through  P  using  equation  M) .  Thus,  S 

is  the  sum  of  the  squares  of  the  differences  between 

Y  ■ 

the  measured  value,  1  ,  and  the  theoretical  value, 
f  (^i  )  at  each  data  point,  ~i  .  We  use  the  vector 

notation  for  ~i  because  there  is  no  necessity  for  the 
data  points  to  be  scalars.  In  fact,  in  this  paper,  the 
ore  the  frequencies  and  the  temperature  and 
y.  is  the  corresponding  measured  imaginary  dielectric  constant 

i. 

In  order  to  use  Newton’s  formula,  we  must  then 

VS 

calculate  ~  and  JR.  The  appropriate  partial 
derivatives  are  calculated  numerically  by  the  following 
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formulae: 


S(a+6a.)  -  S(a-6a. ) 

—  A.JL  »»-l 

26 


and 


92S  (  1  \ 

- fis  (  - -r  I  (S(a+6a-+6a ■)  -  S(a+6a--6a-) 

3a.  3a.  \  4  62  /  ~  1  ~3  -  ^  ^ 

i  3 

-S  (a-6a .  +6a . )  +  S(a-6a  -6a  .) } 

yv  a-  i  ~  -j  /<•  •'*  1 


where  is  a  vector  with  6  in  the  ith  position 

and  zero  elsewhere.  The  starting  values  for  £  are 
supplied  by  the  user.  Also,  the  user  supplies  the 
6~i  ,  which  are  the  sensitivity  values  by  which  the 
parameters  are  varied  in  the  calculation  of  the  partial 
derivatives.  These  sensitivity  values  are  based  on  the 
user's  experience  and  are  often  varied  as  the  fit 
improves.  In  addition  to  this,  parameters  are 
sometimes  " transf orned"  in  order  to  make  them  more 
amenable  to  analysis.  A  case  in  point  is  'o ,  the 
M-ciproca]  frequency  factor.  leather  than  fit  a 
parameter  that  can  easily  vary  by  a  few  orders  of 
nvinitude,  it  is  much  better  to  fit  ln(  o )  and  then 
determine  To  by  exponentiation. 

After  determining  sensitivities,  the  user  then 
commands  the  computer  to  begin  the  least  squares  fit. 


At  the  start  of  the  routine,  the  computer  calculates 
and  stores  S(a)  as  well  as  the  current  parameter 
values.  In  the  process  of  calculating  ,2  s  and  H,  it 

V 

evaluates  S  in  the  vicinity  of  j3,  using  standard 
difference  quotient  approximations  for  the  derivatives. 
If  one  of  these  values  happens  to  be  at  the  current 
minimum  value  for  S,  both  it  and  the  corresponding 
parameter  values  are  stored.  Finally,  ,2s  and  JH  are 
used  to  obtain: 

Aa  =  -  (H)'1 (VS) 

The  values  for 

S(a+Aa)  and  S(a+2Aa) 

are  then  computed.  A  parabolic  fit  is  applied  to  the 
values 

S (a+cAa)  (c  =  0,1,2) 

w  (V 

to  estimate  the  value  c*  for  which 

S  (a+c* Aa ) 

■v  «•' 

is  a  minimum.  Finally, 

S  (a4c*Aa) 

w-> 

is  evaluated.  At  this  point  the  routine  chooses  the 
point,  a*,  for  which  t  tic  computed  value  of  P  was 
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smallest,  and  alters  the  current  parameter  values 
accordingly.  Iterations  are  continued  until  no  further 
significant  progress  is  made. 

There  are  several  valuable  features  to  the 
routine  that  greatly  enhance  the  purely  numerical 
techniques  of  minimization.  The  particular  numerical 
algorithm  chosen  for  doing  the  minimization  provides  a 
nice  balance  between  the  sure,  but  slowly  converging, 
grid  search  method,  and  the  rapidly  converging,  but 
somewhat  unstable  Newton  method.  As  mentioned  in  the 
introduction,  the  parabolic  fit  essentially  recovers 
the  curvature  in  the  Newton-Raphson  (gradient) 
direction  and  thus  describes  to  some  extent  the  metric 
on  the  graph  of  the  squared  difference  sum  function  or 
the  energy  function. 

Finally,  and  by  far  the  most  important  feature 
of  all,  is  the  interactive  nature  of  the  numerical 
techniques  and  graphic  display  of  the  fit  in  its 
current  state.  The  ability  of  the  individual  to 
interpret  and  evaluate  graphical  data  is  far  superior 
to  any  known  purely  numeric  computer  algorithm.  At 
times,  especially  when  the  numerical  methods  fail  to 
improve  the  fit  when  obvious  improvements  can  he  made, 
the  experienced  user  can  manually  adjust  parameter 
values  and  sensitivity  values  in  an  attempt  to  get  the 
parameters  into  a  region  whore  the  numerical  methods 
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will  once  again  converge.  This  interaction  has  proved 
to  be  the  most  valuable  feature  of  all. 


t 


This  interactive  aspect  is  of  course  not 
available  with  the  defect  simulation  package,  but 
experimentation  has  established  that  the  above 
minimization  procedure  nonetheless  leads  to  rapid 
convergence  of  ion  location.  The  instability  where  one 
ion  converges  to  another's  location  (due  to  the  two 
getting  too  close)  can  be  easily  dealt  with  by  giving 
the  ions  an  arbitrarily  large  repulsive  potential 
within  an  appropriate  radius.  Statistics  kept  during 
all  runs  to  date  show  that  both  the  grid  search  and 
parabolic  fit  in  the  Newton  direction  are  techniques 
well  worth  inclusion  in  the  defect  simulation  package. 
Speci f ical ly ,  inclusion  of  the  parabolic  fit  feature 
appears  to  have  sped  up  convergence  by  a  factor  of  ten. 
This  makes  the  running  time  requirements  quite 
acceptable . 

Appendix  C.  Defect  Simulations. 

We  now  give  more  details  of  the  operation  of  the 
defect  simulation  package.  The  model  contains  about 


electrostatic  calculations  for  a  given  ion  in  the 
central  region  are  performed  exactly  for  the 
surrounding  cube  of  llxllxll  fluorines  plus  cations  and 
the  remaining  electrostatic  energy  is  approximated  via 
a  Madelung  sum. 

One  or  more  defect  ions  can  be  added  to  a  file 
containing  the  basic  model  including  substitutionals, 
interstitials,  or  vacancies.  The  lattice  is  then 
allowed  to  relax  to  the  minimum  energy  configuration. 
That  is,  for  each  ion  in  turn  an  approximate  energy 
gradient  is  evaluated,  representing  the  net  force  on 
that  ion  due  to  the  other  ions  in  the  model.  The  ion 
is  allowed  to  move  an  appropriate  amount  in  the 
direction  specified  by  the  force  (gradient),  and  then 
the  computer  moves  on  to  the  next  ion.  All  this  is 
done  a  large  number  of  times  (1000  to  F>oon  moves  being 
necessary)  .  The  energy  function  used  will  be 
discussed  in  the  next  section. 

One  crucial  factor  is  the  question  of  how  far 
each  ion  is  to  move.  This  should  depend  on  the  force 
exerted  on  that  ion  but  should  not  be  too  great,  as  the 
ions  in  actuality  all  move  simultaneously.  One  does 
not  wish  to  teach  a  physically  implausible 
configuration  by  moving  one  ion  around  excessively. 

We  proceed  by  starting  with  initial  values  based 
on  experience  and  trial-and  error  experimentation  with 
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the  model.  Starting  with  a  Newton-Raphson  approach 
loosely  controlled,  reasonable  values  for  the  move  are 
established,  and  these  are  what  are  initially  used. 
Subsequently,  each  ion  is  allowed  to  move  only  as  far 
as  it  moved  the  previous  time;  half  as  far  in  most 
cases  (down  to  a  minimum  size  allowed  move).  Ions  with 
largest  potential  move  are  moved  first.  This  helps 
speed  up  the  convergence  process. 

Having  dealt  with  the  question  of  how  far  to 
move,  the  next  question  is  in  which  direction  each  ion 
is  to  be  moved.  For  each  ion,  the  following  information 
is  stored:  reference  coordinates,  type,  shell 

9 

coordinates,  nucleus  coordinates,  and  move  size.  The 
shell  is  specified  by  a  given  radius  (ion  size)  and  is 
assumed  to  be  roughly  spherical.  Its  coordinates  are 
those  of  its  center,  which  should  be  roughly  the  same 
as  the  nucleus  coordinates.  Move  size  is  the  size  of 
the  largest  distance  the  shell  center  and  nucleus  are 
allowed  to  move. 

The  program  calculates  the  energy 


C'  < x  ]'•••'  x6  > 

stemming  from  the  ion  under  consideration  being  at 


shell  location  (x 


,x  )  and  nucleus  position 


(x,  >x 


) .  The  force  on  the  ion  is  -VC  which  is 


4  5  6 

approximated  using  the  first  difference. 


As  only  the  direction  is  needed,  the  division  by  the 
step  size  6  to  give  the  difference  quotient  is  not 
performed,  to  save  time.  A  standard  Newton-Paphson 
routine  to  locate  a  zero  of  the  force  function  -V G  is 
performed.  Thus  the  force  on  the  ion  is  estimated  to 
be  zero  when  the  ion  at  location  x  is  moved  to  the 

r»* 

new  location  given  by  x+  Ax,  where  Ax  is  given  by 
Ax  =  - (H)-1 (VG) 


Here  H  is  the  Hessian  matrix  of  second  derivatives  of  G 
or  equivalently  the  matrix  of  partials  of  the  force 
F  =  -VG.  See  section  2  above. 

The  second  partials  in  H  are  estimated  as  second 
differences.  Note  that  we  thus  take 

Ax  »  -(H)"l(VG).S 


where  H  is  the  approximate  Hessian  and  VG  the 
approximate  gradient.  In  all  these  difference 
quotients,  division  by  the  step  size  A  can  be 
completely  eliminated,  thus  avoiding  unnecessary 
expensive  machine  divisions.  A  further  benefit  is  the 
increased  numerical  stability  experienced  in  the 
inversion  of  the  matrix  H. 

The  Newton-Raphson  technique  notoriously  suffers 
from  poor  efficiency  when  near  the  zero,  i.e.  when  the 
ion  is  nearly  in  a  minimal  energy  position.  This  is 
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also  true  near  local  minima  or  saddle  points  of  the 
energy  function.  Thus  it  has  proved  useful  to  test  the 
energy  at  the  points  being  used  in  the  partial 
derivative  approximations,  i.e.  conduct  a  "grid 
search",  as  the  energy  has  to  be  calculated  anyway  for 
the  partials. 

Thus  if  the  energy  is  lower  at  one  of  the 
locations  used  in  the  partial  derivative  calculations 
than  at  the  Newton-Raphson  site,  the  former  is  chosen 
as  the  new  ion  location.  A  further  refinement  is 
possible.  It  is  plausible  that  near  the  local  minimum 
of  the  energy  function  (much  as  with  a  least  squares 
function)  that  the  geometry  is  "parabolic",  i.e.  one 
anticipates  that  along  any  line  through  the  current 
location  of  the  ion  the  energy  function's  graph  will  be 
parabolic.  A  parabolic  fit  in  the  Newton-Raphson 
direction  is  performed,  as  described  in  section  ? 
above.  That  most  moves  are  being  made  in  this  fashion 
is  not  too  surprising,  as  the  parabolic  fit  is 
essentially  equivalent  to  finding  the  curvature  at  the 
* ewt on-Raphson  site  in  the  direction  of  the  gradient. 

For  maximum  efficiency,  the  method  of 
approximating  partial  derivatives  is  not  quite  as 
described  above.  First,  note  that  the  matrix  of  second 
partials  of  the  energy  function  (1  for  a  given  ion  may 
be  reasonably  assumed  to  be  symmetric.  Theoretically 
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this  is  justified  as  we  are  working  with  a  matrix  of 
second  partials,  and  continuity  is  not  a  drastic 
assumption. 

The  shell  and  nucleus  in  the  standard  shell 
model  are  joined  by  springs,  i.e.  our  matrix  H  of 
second  partials  is  of  the  form 

i 

-k  0  0  \ 

A  ’  0  -k  n 

0  0  -k  1 

!  i 

0  0 

I  / 

k  0|  B 

0  -k,  / 

Here  k  is  the  spring  constant,  A  is  the  dependence  of 
the  energy  on  the  shell  position,  and  B  is  the  portion 
due  to  the  nucleus. 

We  return  to  the  question  of  approximating  the 
remaining  second  partials  of  the  energy  function.  So 
as  to  combine  a  grid  search  with  a  Newt on-Raphson 
technique,  we  should  estimate  the  partial  ,  for 

example,  as  the  difference  quotient 

( ^  f  2  f  *  *  *  f  ) 

T*\  ~  26 

similarly,  the  second  partial  should  be 

estimated  as 
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32G 


3x1 3x2 


— ^  (G  (x.  +6  ,  x0+6 ,  . . .  ,  x,  )  -G  (x,  +  6 
46  1  ^  b  i 


x2~^ 9  *  *  *  9 ^ ^  ^ 


-G (x^-6 ,x2+6 ,  . . . ,Xg) +G(x1-6 ,x2-6 , . . . ,xg) } 


Initially  <5  differed  for  each  coordinate  hut  it  was 
found  to  be  simpler  to  use  the  same  <$  in  each 
coordinate. 

However,  this  means  that  we  must  evaluate  the 
energy  at  1R  points.  If  one  thinks  of  the  ion  under 
study  as  being  at  the  center  of  a  cube  of  side  ?  , 

we  are  evaluating  at  the  center  of  each  face  and  edge 
of  the  cube.  Instead,  the  program  evaluates  the  enerqy 
at  the  8  points  corresponding  to  the  corners  of  the 
cube  (i.e.  all  of  the  3  shell/nucleus  coordinates 
varied  by  either  +  6  ,  rather  than  just  one  or  two 

being  simultaneously  changed).  The  value  at  the 
middle  of  an  edqe  is  approximated  as  the  average  of  the 
two  corner  values,  for  the  corners  on  the  ends  of  that 
edge.  The  value  of  the  energy  function  at  the  middle 
of  a  face  of  the  cube  is  approximated  as  the  average  of 
the  values  at  the  four  corners  of  that  face.  Thus  at 
the  price  of  a  little  added  complexity,  the  evaluations 
of  the  energy  function  and  hence  the  time  demands  of 
the  program  may  be  halved. 

Note  that  the  previous  move  size  6  is  the  same 
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as  the  _S_  used  in  these  difference  quotients. 

Other  refinements  are:  (i)  the  move  size  6  is 
stored  as  ]>  .001  unless  a  (smaller)  Newton-Raphson  move 
has  bee1"  made.  ( i  i )  if  the  Hessian  is  not  stably 

i 

Invertible,  only  a  "grid  search"  is  performed.  (iii) 

2 

no  division  by  _6_  or  by  6 _  is  performed. 
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